Relaxation Revisited — A Fresh Look at 
Multigrid for Steady Flows 

Thomas W. Roberts*, David Sidilkover^ 
and R. C. Swanson* 


1 Introduction 

The year 1971 saw the publication of one of the landmark papers in compu- 
tational aerodynamics, that of Murman and Cole [9]. As with many seminal 
works, its significance lies not so much in the specific problem that it addressed — 
small disturbance, plane transonic flow — but in the identification of a general 
approach to the solution of a technically important and theoretically difficult 
problem. The key features of Murman and Cole’s work were the use of type- 
dependent differencing to correctly account for the proper domain of dependence 
of a mixed elliptic/hyperbolic equation, and the introduction of line relaxation 
to solve the steady flow equation. All subsequent work in transonic potential 
flows was based on these concepts. Jameson [6] extended Murman and Cole’s 
ideas to the full potential equation with two important contributions. First, 
he introduced the rotated difference stencil, which generalized the Murman and 
Cole type-dependent difference operator to general coordinates. Second, he used 
the interpretation, introduced by Garabedian, of relaxation as an iteration in 
artificial time to construct stable relaxation schemes, generalizing the original 
line relaxation method of Reference [9]. The decade of the 1970s saw an ex- 
plosion of activity in the solution of transonic potential flows, which has been 
summarized in the review article of Caughey [4], 

At about the time of Caughey’s survey, the main thrust of research in com- 
putational aerodynamics was moving away from the full potential equation and 
towards the solution of the steady Euler equations. By analogy with relaxation 
methods for potential equation, solution methods for the Euler equations can 
be thought of as iterations in pseudo-time. Unlike methods for the potential 
equation, by far the most common approach has been to solve for steady flows 
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as the long-time asymptotic solution of the unsteady equations, rather than 
directly solve the steady equations. This is much easier conceptually, as the un- 
steady Euler equations are a hyperbolic system, for which a considerable body 
of theory exists. Abandoning time accuracy allows considerable flexibility in the 
construction of the iterative scheme. To accelerate the convergence to the steady 
state, various types of preconditioning are used. In addition, starting with the 
unsteady equations leads to a straightforward extension of the iterative methods 
to the Reynolds-averaged Navier-Stokes equations, and considerable progress in 
this area has been made in recent years. 

Nevertheless, there is still a great need to improve the convergence rates of 
existing methods. Both the line-relaxation methods for the potential equation 
and the time-iterative methods for the Euler and Navier-Stokes equations suffer 
from slow asymptotic convergence to the steady state. Interpreting the iteration 
as relaxation leads one naturally to consider convergence acceleration methods 
that have been successfully applied to classical relaxation schemes. The foremost 
among such methods is the multigrid algorithm. The theory of multigrid is 
highly developed for elliptic equations, for which 0(n) convergence rates are 
attainable, where n is the number of unknowns in the system. In other words, 
the work required to obtain a solution to the system of equations is proportional 
to the number of unknowns. Classical relaxation schemes for elliptic equations 
are extremely efficient at eliminating the short-wavelength components of the 
error, while the coarse grids in the multigrid process are efficient at removing 
the long- wavelength errors. Application of multigrid acceleration to the Euler 
or Reynolds-averaged Navier-Stokes equations leads one to consider temporal 
integration methods which also provide good damping of the short-wavelength 
error components. 

There are two main classes of multigrid methods based on the unsteady 
equations. One class of methods uses upwind-differencing and implicit time 
integration as the smoother [1, 8, 17]. An alternative approach is one orig- 
inally proposed by Jameson [7]. A finite- volume spatial discretization with 
explicit artificial viscosity is combined with a Runge-Kutta time integration as 
a smoother. This approach has been successfully extended to the Reynolds- 
averaged Navier-Stokes equations [16]. Unfortunately, these approaches have 
resulted in poor multigrid efficiency. When applied to high Reynolds number 
flows over complex geometries, convergence rates are often worse than 0.99 per 
multigrid cycle. Recently, significant improvements have been demonstrated by 
Pierce, et al. [10]. However, when one considers that for the Poisson equation 
on smooth domains convergence rates of nearly 0.1 per cycle are attainable in 
practice, it is clear that there is tremendous room for improvement of existing 
flow solvers. 

In the remainder of this paper, a multigrid algorithm for the Euler equations 
which yields convergence rates comparable to those of the Poisson equation is 
presented. This algorithm abandons the time-marching approach to the steady 
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state, but relies on relaxation of the steady equations. In Section 2, the general 
principles underlying the algorithm are outlined. The mathematical formulation 
of the current approach is given in Section 3, while the solution procedure is de- 
scribed in Section 4. Results for incompressible, inviscid flow in two-dimensional 
channels and around airfoils are shown in Section 5. A brief discussion of the 
extension of the current method to the compressible flow equations is presented 
in Section 6, where a connection with potential flow solvers is also shown. A 
summary is found in Section 7. 


2 An Approach to Multigrid 

According to Brandt [2], one of the major obstacles to achieving ideal multigrid 
performance for advection dominated flows is that the coarse grid provides only 
a fraction of the needed correction for smooth error components. This particu- 
lar obstacle can be removed by designing a solver that effectively distinguishes 
between the elliptic, parabolic, and hyperbolic (advection) factors of the sys- 
tem and treats each one appropriately. The efficiency of such a solver will be 
limited by the efficiency of the solvers for each of the factors of the system. For 
instance, advection can be treated by space marching, while elliptic factors can 
be treated by multigrid. In this example, all components of the error associated 
with the advection terms are eliminated in one sweep, and the convergence rate 
is limited by the speed of the elliptic solver. Brandt presents an approach called 
“distributive relaxation” by which one can construct smoothers that effectively 
distinguish between the different factors of the operator. Using this approach, 
Brandt and Yavneh have demonstrated textbook multigrid convergence rates 
for the incompressible Navier-Stokes equations [3]. Their results are for a sim- 
ple geometry and a Cartesian grid, using a staggered-grid discretization of the 
equations. 

In a closely related approach, Ta’asan [15] presents a fast multigrid solver 
for the compressible Euler equations. This method is based on a set of “canon- 
ical variables” which express the steady Euler equations in terms of an elliptic 
and a hyperbolic partition [14]. Ta’asan uses this partition to guide the dis- 
cretization of the equations. A staggered grid is used, with different variables 
residing at cell, vertex, and edge centers. In Reference [15] it is shown that 
ideal multigrid efficiency can be achieved for the compressible Euler equations 
for two-dimensional subsonic flow using body-fitted grids. One possible lim- 
itation of the use of canonical variables is that the partition of the inviscid 
equations is not directly applicable to the viscous equations. 

Recently the authors [11] have presented an alternative scheme to Brandt’s 
distributive relaxation and to Ta’asan’s canonical variable decomposition. This 
scheme does not require staggered grids, but uses conventional vertex-based 
hnite- volume or Unite-difference discretizations of the primitive variables. This 
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simplifies the restriction and prolongation operations, because the same opera- 
tor can be used for all variables. A projection operator is applied to the system 
of equations, resulting in a Poisson equation for the pressure. By applying 
the projection operator to the discrete equations rather than to the differential 
equations, the proper boundary condition on the pressure is satisfied directly. 
The Poisson equation for the pressure may be treated by Gauss-Seidel relax- 
ation, while the advection terms of the momentum equation are treated by 
space-marching. Because the elliptic and advection parts of the system are de- 
coupled, ideal multigrid efficiency can be achieved. Compared to distributive 
relaxation and the canonical variables approaches, this method is extremely 
simple. 


3 Mathematical Formulation 

The incompressible Euler equations in primitive variables are 

uu x + vu y + p x = 0, 

UV X + VVy + Py = 0 , 

U X + Vy = 0, 

where u and v are the components of the velocity in the x and y directions, 
respectively, and p is the pressure. The density is taken to be one. The advection 
operator is defined by 


Q = ud x + vd y , (1) 

where d x , d y are the partial differentiation operators. The Euler equations may 
be written as 


(Q 0 d x \ fu\ 

Lq = 0 Q d y u = 0. (2) 

\d x dy 0 J \pj 

Introducing the adjoint to Q, defined by 

Q*(f) = —d x (uf) - d y (vf), (3) 

a projection operator P is defined: 

// 0 0 \ 

P = 0 / 0 . (4) 

\d x dy Q*J 
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Applying the projection operator to the Euler equations yields 

(Q 0 d x \ (u\ ( o \ 

Lq = PLq = 0 Q d y U + 2 0 , (5) 

\0 0 A/ \pj \(d x v)(dyu) - (d x u)(d y v) ) 

where A is the Laplacian. The matrix operator on the right-hand side consists 
of the principal part of L (i.e., the highest-order terms of the operator), and 
the remaining terms are are the subprincipal terms. These terms arise because 
the coefficients u and v in the operators Q and Q* are not constant. It is 
important to note that the subprincipal terms can be ignored for the purpose 
of constructing a relaxation scheme. 

The system of equations (5) is a higher-order system than the original Euler 
equations (2). The continuity equation, which is a first-order partial differential 
equation, has been replaced by a second-order differential equation for the pres- 
sure. One might expect that Eq. (5) would require a boundary condition on the 
pressure in addition to the physical boundary condition of flow tangency at the 
wall, which is required by Eq. (2). However, at the boundary of the domain, 
the third equation of (5) takes the form 

(- udyV + vdyU + d x p)h x + (ud x v - vd x u + d y p)h y = 0, (6) 

where h X} h y are the components of the unit normal at the wall. This is simply 
the equation for the momentum normal to the wall. Because the pressure equa- 
tion at the wall takes the form of Eq. (6), which in this case may be thought of 
as a compatibility condition of the governing equations, no auxiliary boundary 
condition on the pressure is needed. 

The operator on the left-hand side of Eq. (5) is upper triangular. Because 
the pressure satisfies a Poisson equation a conventional relaxation method, such 
as Gauss-Seidel, can be used to solve it. Upwind differencing of the advection 
operator in the momentum equations and a downstream ordering of the grid 
vertices allows marching of the momentum equations. A collective Gauss-Seidel 
approach is used here, where the vertices are ordered in the flow direction. This 
is described more fully in the next section. 

4 Solution Procedure 

The first step in approximating L is to discretize the Euler equations (2). Un- 
like the methods of References [3, 15], which use staggered grids, the current 
approach is vertex-based, where all the unknowns are stored at the vertices of 
the grid. Discretizations for quadrilateral structured grids and triangular un- 
structured grids have been coded. A great deal of flexibility in the form of the 
discrete approximation to the momentum equations is possible with the current 
method. 
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1 ► x 

Figure 1: Triangular cell of an unstructured grid 


By way of illustration, consider a cell-vertex discretization on an unstruc- 
tured triangular grid. A typical grid cell 0 is shown in Fig. 1. Cell-averaged 
gradients of the unknowns are found by using a trapezoidal rule integration 
around the boundary of 0. For example, the discrete approximation to the 
gradient of u on 0 is 

iC " ' }<>'." = :r4- {MVi ~ 2 / 2 ) + UiiV'2 ~ yo) + MVo ~ 2 / 1 )) 

“7 ( 7 ) 

- r~r~ («o(-?’l - ^2) + - -?’o) + U 2 (x 0 ~ 

where Aq, is the area of the triangle. The. superscript h is used to denote the 
difference approximation to the corresponding differential operator. Gradients 
of v and p are obtained likewise. These gradients are used to approximate 
Ecp (2) on the triangle. An upwind approximation to Q at the vertices of the 
grid is obtained by distributing the cell-averaged momentum equation residuals 
to the vertices of each triangle appropriately. The current scheme is not tied to 
any particular form of the upwind discretization. One choice, which was used 
to obtain the unstructured grid results presented in Section 5, is the advection 
scheme of Giles, et ah [5]. In this scheme, the residuals are distributed to the 
vertices of 0 using the weights 



where (' n is the length of the projection of 0 onto the crossflow direction, and A n t 
is the component of the length in the crossflow direction of the edge opposite 
the z-th vertex. An alternative upwind discretization currently being devel- 
oped by the authors is based on the multidimensional upwind formulation of 
Sidilkover [12]. 

Once the cell-averaged residuals of the continuity and momentum equations 
have been computed, the projection operator P is applied to these discrete 
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equations to obtain the residual for the pressure Poisson equation of Eq. (5). 
Letting R Pi be the pressure equation residual at vertex i, the application of P 
can be written in integral form, 



where Ai is the area of the control volume centered on the i-th vertex and the 
superscript h is used to denote the discrete approximations to the corresponding 
differential operators on the triangle as before. The summation is over all the 
triangles adjoining the i-th vertex. The dashed lines in Fig. 1 are the segments of 
the boundaries of A 0 , Ai and A 2 that lie in cell 0. The boxed terms in (8) are the 
cell-averaged residuals of the x and y momentum equations and the continuity 
equation. The contributions of the cell-averaged residuals on 0 to R Pi are found 
by evaluating Eq. (8) over the appropriate segment of the boundary of Ai lying 
in 0, taking the boxed terms in Eq. (8) to be constant over the cell. 

Applying the projection operator P at the discrete level in this way, rather 
than starting with the differential equations (5) and discretizing them, has two 
important advantages. First, the discrete approximation of Eq. (8) at bound- 
ary vertices reduces to Eq. (6), automatically providing the correct boundary 
condition for the pressure. Second, if the momentum and continuity equations 
are discretized on the triangles in conservation form, it is possible to obtain a 
fully conservative scheme. This is particularly important for compressible flows 
with shocks. 

The multigrid algorithm uses a sequence of grids Gk, Gk- i, . . . , G o, where 
Gk is the finest grid and Go the coarsest. Call the discrete approximation to 
the operator L on the &-th grid L*,, and let be the solution on that grid. 
This system has the form = f^, where the entries of are 3x3 block 

matrices which operate on the unknowns (it, v,p) T at each grid vertex. A general 
iteration scheme is constructed by writing the operator as — N*,, 

where the splitting is chosen such that is easily inverted. Lexicographic 
Gauss-Seidel is obtained by taking to be the block lower-triangular matrix 
resulting from ignoring the terms above the diagonal blocks of L*,. A further 
simplification is obtained if the diagonal blocks of contain only those entries 
corresponding to the principal part of the operator. Because the operator in 
Eq. (5) is upper triangular the diagonal blocks of will then be 3 X 3 upper 
triangular matrices. 

Letting qjj be the re-th iterate of the solution on the £;-th grid, the relaxation 
iteration is 

M fc q £ +1 = f k + N k q n k . 


7 



The operator is nonlinear, so 1VR and N*, are functions of q k and q k +l ■ 
Letting Sq k = q k +l — q]J, the iteration may be rewritten as 

M fc hq£ = f k - l k q n k . (9) 

Because 1VR is block lower-triangular, hqJJ is found by forward substitution. At 
each vertex, a 3 X 3 upper triangular matrix must be inverted. 

If the discrete approximation to the advection operator Q is fully-upwind 
and the grid points are ordered in the flow direction, then the 3x3 blocks of 
will have zeroes in the first two rows. In this case, lexicographic Gauss-Seidel 
relaxation is equivalent to space-marching of the advection terms. The advected 
error is effectively eliminated in one relaxation sweep and the convergence rate of 
the system becomes that of the Poisson equation for the pressure. It is possible 
to get ideal multigrid convergence rates because each component of the error is 
treated appropriately. 

A straightforward Full Approximation Scheme (FAS) multigrid iteration is 
applied to the system of equations. Let ~L k -i be the coarse grid operator, I k _ x 
be the hne-to-coarse grid restriction operator, and I k ~ x be the coarse-to-hne 
grid prolongation operator. If cp, is the current solution on grid £;, the residual 
on this grid is r k = i k — L^cp,. This leads to the coarse-grid equation 

Lfc-iqfc-i = ffc-i = I k _i rfc + Lfc-i {l k -\ qfc) • (10) 

After solving the coarse-grid equation for q^-i, the fine-grid solution is corrected 
by 


qr <- q* + It 1 (q*-i - iLA) • (n) 

Equation (10) is solved by applying the same relaxation procedure that is used 
to solve the fine-grid equation. Multigrid is applied recursively to the coarse-grid 
equation. On the coarsest grid, many relaxation sweeps are performed to insure 
that the equation is solved completely. A conventional R-cycle or W - cycle is 
used. 


5 Results 

Both unstructured grid and structured grid flow solvers based on the theory in 
Sections 3 and 4 have been written. These codes are described in Reference [11], 
where extensive solutions are presented. Results illustrating the efficiency of the 
scheme are presented here. 

Solutions for incompressible, inviscid flow in a channel have been obtained 
with both solvers. The channel geometry and boundary conditions are shown 
in Fig. 2. The shape of the lower wall between 0 < x < 1 is y(x) = rsin 2 7rx. 
For the computations shown here, the thickness ratio r is 0.05. The flow angle 
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Figure 2: Channel geometry. 


and total pressure are specified at the inlet and the pressure is specified at the 
outlet. The flow tangenc.y condition u-n = 0 is enforced at the upper and lower 
walls of the channel. Solutions were obtained on quasi-uniform quadrilateral 
grids. A simple shearing transformation was used in the center part of the 
channel to obtain boundary conforming grids. For the unstructured grid solver, 
the grids were triangulated by dividing each quadrilateral cell along a diagonal. 
A series of nested coarse grids was obtained by coarsening the fine grids by a 
factor of two in each coordinate direction. In all cases shown below, the coarsest 
grid was 7x3 vertices. Lexicographic Gauss-Seidel relaxation was used, with 
the grid vertices ordered from the lower-left to the upper-right of the channel. 
This resulted in downstream relaxation of the momentum equations. A V(2, 1) 
multigrid cycle was used; that is, two relaxation sweeps were performed on 
each grid before restricting to the coarse grid, and one relaxation sweep was 
performed after the coarse-grid correction was added to the fine-grid solution. 

The computed pressure on a grid of 97 X 33 vertices is shown in Fig. 3 
for the unstructured grid flow solver and in Fig. 4 for the structured grid 
solver. Comparisons of convergence rates for different grid densities are shown 
in Figs. 5 and 6 for the unstructured and structured grid flow solvers, respec- 
tively. Thq Xi norm of the pressure equation residual is shown; the momentum 
equation residuals show the same behavior. The finest grid used for each flow 
solver contained 385 X 129 vertices, with a total of 7 grid levels. 

The convergence rate of the unstructured grid solver on the finest grid is ap- 
proximately 0.190 residual reduction per multigrid cycle. The structured grid 
results are slightly better at 0.167 per cycle. These rates are comparable to the 
ideal rate of 0.125 per cycle for the Poisson equation. The better performance 
of the structured grid solver is most likely because of better restriction and 
prolongation operators; the unstructured flow solver performs bilinear interpo- 
lation using only the locations of a fine-grid vertex and the three vertices of 
the coarse-grid cell containing that vertex. What is most important is that the 
figures show nearly ideal multigrid convergence rates, independent of the grid 
spacing. This shows that convergence is achieved in order n operations. 

For complex geometries it may not be practical to generate a series of nested 
unstructured grids, and the performance of the multigrid solver may be expected 
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Figure 3: Pressure, contour increment A p = 0.01, for an unstructured grid 
of 97 x 33 vertices. 



Figure 4: Pressure, contour increment A p = 0.01, for a structured grid of 97 x33 
vertices. 

to deteriorate. To show the robustness of the current method, the triangular 
grid solver was run for a series of non-nested coarse grids. These were generated 
by randomly perturbing the locations of the vertices on each of the nested grids 
independently. The perturbed 49 x 17 grid is shown in Fig. 7. The computed 
pressure on a perturbed 97 x 33 fine grid with 5 grid levels is shown in Fig. 8 and 
the convergence rate is shown in Fig. 9. The pressure contours are very smooth, 
showing no sign of the lack of grid smoothness. The asymptotic convergence 
rate has deteriorated to a still-respectable 0.24 per cycle. 

Solutions for nonlifting flow over a symmetric Karman-Trefftz airfoil have 
been obtained with the structured grid solver. A fine O-grid of 385 X 193 vertices 
was generated from a conformal mapping, and the coarse grids are nested by 
recursively eliminating every other vertex in each coordinate direction. The grid 
spacing was chosen to obtain unit aspect ratio grid cells. The outer boundary is 
approximately 13 chord lengths from the airfoil. Far-held boundary conditions 
are given by the analytic solution. At inflow points along the outer boundary 
the total pressure and flow inclination angle are specified. For outflow points the 
pressure is specified. On the airfoil surface the tangency condition is enforced. 

To obtain ideal multigrid convergence rates, it is necessary to sort the ver- 
tices in a downstream order so that the advection terms in the momentum 
equations are marched. This is easily done here by relaxing along the radial 
grid lines from the outer boundary to the airfoil surface over the forward half of 
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V(2, 1) cycles 

Figure 5: Comparison of convergence rates on unstructured grids 



V"(2,l) cycles 


Figure 6: Comparison of convergence rates on structured grids. 






Figure 7: Grid generated by perturbing the vertices of the 49 x 17 grid. 



Figure 8: Pressure, contour increment A p = 0.01, randomly perturbed unstruc- 
tured grid of 97 x 33 vertices. 


the domain, and from the airfoil surface to the outer boundary over the latter 
half of the domain. For each case run, the coarsest grid consisted of 13 X 7 
vertices. 

Comparisons between computed and analytic surface pressure coefficients 
for nonlifting flow around the Karman-Trefftz airfoil are shown in Fig. 10. 
A 1 / F(2,1) multigrid cycle was used for these computations. The computed 
solution agrees very well with the analytic solution, except for the recompres- 
sion at the trailing edge. Note that there is no clustering of the grid in this 
region, which exacerbates the problem. 

A comparison of the convergence rates of the pressure equation residual for 
three grid densities is shown in Fig. 11. A slight deterioration of the convergence 
rate with increasing grid refinement is observed: on the 385 X 193 grid, the 
rate is 0.153 per cycle. Nevertheless, as with the channel flow results, the 
convergence rates are very nearly grid independent, and are very close to the 
ideal rate of 0.125 per cycle. 

A summary of the convergence rate on the finest grids is presented in Table 1. 
Two sets of results are shown: the convergence rate per multigrid cycle, and 
the convergence rate per work unit. For the purposes of the discussion a work 
unit (WU) is taken to be one Gauss-Seidel relaxation sweep on the finest grid. 
This is essentially the cost of one residual evaluation on the finest grid. The 
actual convergence rates are compared to the ideal convergence rates, which are 
computed as follows. Let fi be the smoothing rate of the relaxation method. For 
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Figure 9: Convergence rate, randomly perturbed unstructured grid of 97 X 33 
vertices, 5 grid levels, V(2, 1) cycle. 


a V(m,n) or lU(m,n) cycle, the ideal convergence rate is fj, m+n . Lexicographic 
Gauss-Seidel for the Poisson equation has a smoothing rate fi = 0.5. This gives 
an ideal convergence rate of 0.5 3 = 0.125 for a (2,1) or a W{ 2,1) cycle. To 
compute the convergence rate per work unit we use the following formula. For 
each cycle, there are a total of m + n fine-grid relaxation sweeps. Examination 
of Ecp 10 shows that the hne-to-coarse grid restriction requires one residual 
evaluation on the fine grid and an additional residual evaluation on the coarse 
grid. Thf coarse grid residual evaluation is 1/4 the cost of a fine grid residual 
evaluation. Because most of the cost of a relaxation sweep is in the evaluation 
of the residual, we have that each cycle requires a total of (m + n + 1 + 1/4) work 
units on the finest grid. The cost of interpolating the residuals and solutions 
between grid levels is neglected. 

For a V(m, n)-cycle, we have that 


WU 1 1 W1 1 1 

__ K(m + , ! + l + _)( 1+ _ + _ + ...) 

4, 5 

3 4 ; 

Because a IT-cycle involves two coarse-grid solutions per cycle, we have 

WU 1 1 W1 1 1 

— — ss m + n + 1 + - 1 + - + 7 4 

IT-cycle 4 2 4 


- 'J , 

= 2 ( 7 ??. + n -| — ). 

4 


These numbers yield ideal convergences rates of /t< 3 * m+n )/( 4 ( m + n + 5 / 4 )) p er work 
unit for a V-cycle and /U ( m + n )/( 2 ( m + n + 5 / 4 )) p er work unit for a IT-cycle. 
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X 


Figure 10: Surface pressure coefficient, nonlifting Karman-Trefftz airfoil, 193 X 
97 grid. 



Figure 11: Comparison of convergence rates for nonlifting Karman-Trefftz air- 
foil. 
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Convergence Rate 

Case 

Cycle 

per 

cycle 

per work unit 



ideal 

actual 

ideal 

actual 

channel, 385 X 129, unstructured 

V(2, 1) 

0.125 

0.190 

0.693 

0.746 

channel, 385 X 129, structured 

V'(2. 1) 

0.125 

0.167 

0.693 

0.729 

airfoil, 385 X 193, structured 

W(2,l) 

0.125 

0.153 

0.783 

0.802 


Table 1: Summary of convergence rates for multigrid solver on finest grids for 
channel and airfoil flows, with a comparison to the ideal rates. 


The V(2, 1) cycle is seen to require 5% WU per cycle. The W(2, 1) is 50% 
more expensive, requiring WU per cycle. By way of comparison, one V(2, 1) 
cycle is only slightly more work than a single time step of a 5-stage Runge-Kutta 
scheme on the finest grid. The ideal convergence rates for lexicographic Gauss- 
Seidel is of 0.693 per WU for a U(2, 1) cycle and 0.783 per WU for a W(2, 1) 
cycle. The actual rates shown in Table 1 are seen to be very close to ideal. 

The convergence rates in Table 1 can be used to estimate the work required 
to obtain a solution to the level of the discretization error on the fine grid. Let p 
be the order of approximation of the discrete operator and let hk be the grid 
spacing parameter on the Uth grid. An initial guess to the solution on the fine 
grid Gk is obtained by interpolating a solution computed on grid Gk- i- Assume 
that the solution on Gk- i has been obtained to the level of the discretization 
error t r _ 1 = 0(h p ) on that grid. The multigrid cycle is used to reduce the 
error from t r to t r . Letting /% be the convergence rate per work unit, the 
amount of work Wk required to get the solution on Gk from the initial solution 
on Gk- i is 


W K = 


i°g/% 


log 


P 


‘ K - 1 


lo g/% 


log 


h 


K—l 


A Full Multigrid (FMG) cycle starts with a solution on the coarsest grid, Go, and 
recursively generates improved solutions on the finer grids using the strategy 
above. 

For the nested grids considered here, the grid spacing parameters are related 

by hk - 1 = 2 hk, and the amount of work on each grid is related by Wk - 1 = 

Wk/ 4. The discretization is second-order accurate, i.e., p = 2. This gives us the 

estimate for the total work to obtain a solution accurate to r„ to be 

A 

U'total = j log ( A 

lo g/V V 2 / 

8 log 2 

3 log/% 


1 + j + ii + "' 


( 12 ) 
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Using the values in the last column of Table 1 for p W7 we see that channel flow 
solutions can be obtained to the level of discretization error in approximately 
6.4 WU using a FMG cycle. Airfoil solutions can be obtained in about 8.3 WU. 
These estimates are generally low, and in fact are less than work of a single FMG 
cycle (7.6 and 11.3 WU for the channel and airfoil cases, respectively). The work 
computed using Eq. (12) also does not account for the introduction of short- 
wavelength errors in the interpolation of the coarse grid solutions to the fine 
grids. Nevertheless, Eq. (12) is a useful guide to the expected performance of 
the multigrid scheme. 


6 Extension to Compressible Flow 


The scheme presented here has a straightforward extension to the compressible 
Euler equation. In primitive variables the equations are 


Lq = 


o 

o 

o 


M 

0 pQ 0 d x 


u 

0 0 pQ d y 


V 

\0 pd x pd y jzQJ 


w 


= 0, 


where p is the density, c is the speed of sound, and s is the entropy, 
projection operator P for this system is 


(l 

0 

0 

Vo 


0 

I 

0 

d T 


0 

0 

I 

d v 


o\ 
0 
0 

Q'J 


(13) 


The 


(14) 


The operators Q and Q* are defined as in Eqs. (1) and (3). Applying this to 
Eq. (13) and ignoring the subprincipal terms as before yields 


PLq = 
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0 pQ 0 d x 


u 

0 0 pQ d y 


V 

\0 0 0 A — M 2 d 2 s j 


w 


+ s.p.t., 


(15) 


where M is the Mach number, d s is the partial derivative in the streamwise 
direction, and “s.p.t.” are the subprincipal terms. 

The most significant difference between the compressible and the incom- 
pressible equations is that a Prandtl-Glauert-like operator acts on the pressure. 
Note that this system approaches the system for the incompressible equations 
in the limit of vanishing Mach number. For subsonic flow the compressible 
equations can be solved by the same relaxation scheme as the incompressible 
equations. Unlike time marching methods, the convergence rate will not dete- 
riorate as the Mach number approaches zero. 
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The appearance of the Prandtl-Glauert operator in the pressure equation is 
significant. In effect, the problem of solving the pressure equation is no different 
than that of solving the full potential equation. This is a two-edged sword. On 
the one hand, the difficulties of relaxing the pressure equation in the transonic 
case are precisely those of relaxing the transonic full potential equation. This 
is a fundamental difficulty which is faced by any method that works directly on 
the steady flow equation. On the other hand, one can expect that the wealth 
of experience in solving potential flows can be directly applied to the current 
scheme for the compressible Euler equations. The treatment of the advection 
terms is not essentially different from the incompressible case. 

7 Conclusions 

Murman and Cole introduced type-dependent differencing and relaxation meth- 
ods into computational aerodynamics; practical and efficient methods for solving 
nonlinear flow equations were the result. In subsequent years, the emphasis has 
shifted toward iterative methods based on the unsteady equations. In this pa- 
per, it has been shown that great improvements in the efficiency of flow solvers 
can be achieved by changing the point of view from the unsteady to the steady 
equations. As Murman and Cole introduced type-dependent differencing, so the 
current method relies on a discretization which distinguishes between the ellip- 
tic and hyperbolic parts of the system. As Murman and Cole used relaxation to 
solve the steady equations, so the current method applies relaxation with multi- 
grid to the steady equations. This approach yields textbook multigrid efficiency 
for the steady Euler equations. It is a particularly simple approach; conventional 
finite-difference or finite-volume discretizations of the governing equations may 
be used, allowing flexibility in the choice of the underlying numerical method. 
Unlike time-marching approaches, but like potential flow methods, the conver- 
gence rate of the method does not degrade for low-speed flows, and the correct 
incompressible limit is recovered. Finally, this method can be applied to incom- 
pressible, viscous flow following the ideas of Sidilkover and Ascher [13]. 

There remains a great deal of work to be done before the full Reynolds- 
averaged Navier-Stokes equations can be solved as efficiently as the simple prob- 
lems shown here. The present work only addresses one particular, but never- 
theless important, aspect of the problem, namely the appropriate discretization 
of the governing equations. If textbook multigrid efficiency is to be achieved for 
compressible, viscous flow, it will likely require an approach along the general 
outlines presented here. 
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